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Abstract — Methods to recover jump-sparse and sparse signals 
from blurred Incomplete data corrupted by (possibly non- 
Gaussian) noise are developed. These reconstruction tasks are 
approached by the minimization of Inverse Potts functionals. 
Analytical results (existence of mlnlmlzers, complexity) on inverse 
Potts problems are obtained and relations to sparslty problems 
are provided. A new method for minimizing Inverse Potts 
problems, which Is based on dynamic programming and the 
alternating direction method of multipliers (ADMM), Is proposed. 
It Is shown In a series of experiments that the proposed method 
yields very satisfactory jump-sparse and sparse reconstructions, 
respectively. The capability of the method Is highlighted by 
comparing It with classical and recent approaches such as 
TV minimization (jump-sparse signals), orthogonal matching 
pursuit, Iterative hard thresholding, and iteratively rewelghted 
£^ minimization (sparse signals). 

Index Terms — Jump-sparslty, sparslty. Inverse Potts functional, 
piecewise constant signal, denolslng, deconvolutlon, incomplete 
data, segmentation, ADMM. 

I. Introduction 

In this article we aim at reconstructing jump-sparse (and 
sparse) signals x £ M" from linear noisy measurements b E 
M™ (or C") given by 

b = Ax + noise, 

where A is a (general) m x n matrix. The reader may 
think of A being a Toeplitz matrix modeling the blur or a 
Fourier matrix, or a combination of both. In particular, we 
deal with possibly incomplete data meaning that the number 
of measurements m is significantly smaller than the size 
n of the original signal. Jump-sparse signals, i.e., signals 
which only jump at a few not a priori known locations, arise 
in various applications. Examples are stepping rotations of 
bacterial flagella [1], the cross-hybridization of DNA [2], [3], 
[4], single-molecule fluorescence resonance energy transfer 
[5], and MALDI imaging [6]. 

The above reconstruction problem is in general ill-posed 
and therefore needs regularization. This is usually achieved 
by minimizing a suitable energy functional which expresses 
a tradeoff between data-fidelity and regularity (Tikhonov-type 
regularization). In view of the jump-sparsity of the underlying 
signal, the number of jumps ||Va;||o is a natural and powerful 
regularizing term. The corresponding minimization problem, 
called inverse Potts problem (iPotts), reads 



P^{x)^j-\\Vx\\o + \\Ax-b\\P 



(1) 
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Here the parameter 7 > controls the tradeoff between jump- 
sparsity and data fidelity, which is measured by some £p norm, 
p > 1. If the noise is Gaussian then p = 2 is the natural choice 
whereas our experiments indicate that for more general noise, 
such as Laplacian or impulsive noise, p = 1 is the better 
choice. 

The inverse Potts functional is not convex. To avoid the 
resultant difficulties, the jump-penalty |lVx||o is often replaced 
by the convex total variation ||Va;||i = X^i l^«+i ^ ^i\- 
This substitution guarantees access to the tools of convex 
optimization. Fast algorithms (e.g., [7], [8], [9]) are available 
and the theory guarantees convergence to a global minimum. 
However, the minimizers of the TV problem in general differ 
from those of the inverse Potts problem; see for example 
Figure 1 . It turns out that minimization of the Potts functional 
yields genuine jump-sparse signals whereas TV minimization 
does so only approximately. 

In this work, we are first concerned with the question of 
existence of minimizers which is more involved than it seems 
at first glance. In fact, we will see that the finite dimensional 
inverse Potts problem (1) has a minimizer whereas its contin- 
uous time counterpart in general need not have a minimizer 
at all. We further show that the inverse Potts problem is NP- 
hard; thus we cannot compute exact minimizers efficiently. 
Accepting this fact, we develop a strategy based on ADMM 
which shows good performance in practice. Furthermore, we 
shed light on the relation between the jump-sparsity problem 
(1) and the sparse recovery problem. Let us be more precise. 

A. Proposed ADMM approach to the inverse Potts problem 

We approach the inverse Potts problem (1) using the al- 
ternating direction method of multipliers (ADMM). ADMM 
algorithms were successfully applied to convex optimization 
problems [10], [11], and, in particular, to TV minimization 
[12], [13]. We propose the iteration 
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where the parameter /i^, is updated by ^k+i = t ■ ^k with 
fixed T > 1. The key point is that each subproblem of 
this iPotts-ADMM algorithm is numerically tractable. The first 
one is a classical Potts problem (equation (1) with A = id) 
which can be solved fast and exactly in the univariate case. 
For multivariate data such as images we use the satisfactory 
heuristics of [14], [15]. The second subproblem consists of 
minimizing a classical Tikhonov functional. When p = 2, we 



solve a normal equation and, for p = 1, we use a fast semis- 
mooth Newton method [16]. As the inverse Potts problem is 
NP-hard, we cannot expect that the method provides global 
minimizers of (1) for all input data. Nevertheless, we show 
that our algorithm converges. In the numerical experiments, 
we obtain very satisfactory reconstruction results. 

B. Inverse Potts problems and sparsity 

The inverse Potts problem is closely connected to the 
"Lagrangian formulation" of the sparse recovery problem 



Sj{x) 



-f ■ \\x\\o + \\Ax - b\\P 



(2) 



The formulation (2) has been considered in [17], [18], for 
instance. General references concerning sparsity are the books 
[19], [20], [21] where also a variety of applications may be 
found in. 

As with the inverse Potts problem and TV minimization, 
one can replace the number of jumps ||a:||o by the absolute 
sum X^il^il '^o obtain a convex relaxation of the sparsity 
problem (2) called basis pursuit denoising (BPDN) or fi- 
minimization. It is one topic of compressed sensing [22], 
[23], [24] to clarify under which conditions a minimizer of 
the €^ -functional minimizes the sparsity problem (2). Positive 
answers (with a high probability) are obtained under quite 
restrictive assumptions on the matrix A such as the restricted 
isometry property [25]. If such conditions are not met the 
solutions of BPDN are in general not minimizers of (2). 
Further related work replaces the jump-penalty ||a;||o by the 
non-convex functionals ||a;||^ with < q < 1 [18], [26], [27]. 

In this work, instead of using relaxations, we reduce the 
sparsity problem (2) to an inverse Potts problem of the form 
(1). We show that this can be done for all data fidelity terms 
based on the p-norm with p > 1. Thus we may approach 
the sparsity problem (2) using the proposed iPotts-ADMM 
algorithm. We further obtain that the inverse Potts problem 
can be cast into a sparsity problem for euclidean data terms 
(P = 2). 

C. Applications and numerical experiments 

We apply the proposed iPotts-ADMM algorithm to re- 
construct signals having few jumps from indirect and noisy 
measurements. More concretely, the indirect measurements in 
our experiments are either Fourier data or blurred data which 
both may be incomplete. Thus the measurement matrix A is 
either a reduced Fourier matrix or a reduced Toeplitz matrix. 
The noise in our examples is either Gaussian noise, Laplacian 
noise, or impulsive noise. We observe that the iPotts-ADMM 
algorithm is capable of recovering jump sparse signals almost 
perfectly, for a reasonable level of noise. We compare the 
results of our iPotts-ADMM algorithm with those obtained by 
total variation minimization. We further compare the inverse 
Potts functional and the classical Potts functional (equation 
(1) with A = id) with regard to the segmentation of blurred 
images. 

In addition, we apply the iPotts-ADMM based method to 
sparse recovery problems. As for jump sparse signals we 



consider blurred data under different types of noise. In our nu- 
merical experiments, we achieve similarly good results as for 
jump-sparse signals. We highlight the capability of our method 
by comparing it with orthogonal matching pursuit [28], [29], 
[30], basis pursuit denoising [31], iterative hard thresholding 
[17] and iteratively reweighted £^ minimization [32], which 
are the state-of-the-art approaches to sparse recovery. 

In order to guarantee reproduciblity an implementation of 
our algorithms is freely available at http://pottslab.de. 

D. Outline of the paper 

We start out to formulate our theoretical results on the 
inverse Potts problem in Section II. In Section III, we derive an 
ADMM algorithm for the inverse Potts problem. In Section IV 
and Section V, we provide numerical experiments; Section IV 
deals with jump-sparse signals whereas, in Section V, we 
consider sparse signals. Finally, we supply the proofs in 
Section VI. 

II. Inverse Potts problems and their relation to 

SPARSITY 

We start our analysis of the inverse Potts problem by 
considering the question of existence of minimizers. It is 
remarkable that there is a significant difference between the fi- 
nite dimensional discrete time case and its infinite dimensional 
continuous time counterpart. More precisely, we obtain a 
positive answer for the discrete time problem (1) but a negative 
answer for the corresponding continuous time problem. 

Theorem 1. The inverse Potts problem (1) has a minimizer. 

The proof of Theorem 1 is given in Subsection VI-A. It 
uses the compactness of the closed unit ball and the lower 
boundedness of an injective linear mapping which are features 
of finite dimensional spaces. Thus it does not carry over to 
the infinite dimensional continuous time case. Moreover, the 
next theorem states that the continuous time counterpart of 
Theorem 1 is false in general. The continuous time counterpart 
of (1) is obtained by replacing the finite dimensional signal 
and data spaces by L^ function spaces and the matrix A by a 
bounded operator A between those function spaces. 

Theorem 2. There are linear operators A and data b in 
LP such that the continuous time inverse Potts problem with 
respect to A and b does not have a minimizer 

The proof of Theorem 2 is given in Subsection VI-A. The 
explicit counter-examples we give are convolution operators 
which are in fact important from a practical point of view. 

The next natural step after showing the existence of min- 
imizers (in the dicrete case) is to clarify the complexity of 
computing such a minimizer We obtain the following result. 

Theorem 3. The inverse Potts problem (1) is NP hard. 

As a consequence, a fast exact algorithm is not available 
(unless P — NP) and one has to resort to approximative 
strategies (see Section III). The proof of Theorem 3 is given 
in Subsection VI-B. 
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Fig. 1: Tlie original signal (dashed line) is convolved by a Gaussian kernel of standard deviation 6 which m = 128 measurements are randomly selected from. 
The resulting data is corrupted by different types of noise (left). The total variation method (TV) mainly reconstructs the constant parts but adds transitional 
points in between the plateaus for Gaussian and Laplacian noise. The iPotts-ADMM recovers the true signal almost perfectly; in particular, the correct number 
of jumps. For impulsive noise, the iPotts-ADMM and the TV method perform equally well. 



Finally, we are interested in the relations between sparsity 
problems and inverse Potts problems. We start with the easy 
direction and first consider the sparsity problem (2). We 
find a corresponding inverse Potts problem whose minimizers 
are directly related to the minimizers of the initial sparsity 
problem. 

Theorem 4. Let x* € M"+^ be a minimizer of the inverse 
Potts functional associated with the matrix B = AV , i.e.. 



B, b',i. e., 



X* £ arg min 7||Va;||o 



\Bx 



%■ 



(3) 



Then u* = Va;* minimizes the sparsity problem (2) related to 
the matrix A and data b. 

We obtain a converse result for p = 2. The relations between 
the matrices A and B and between the data can still be given 
explicitly but are not as simple as above. Our construction 
does not work for general p 7^ 2 and it is not clear to us how 
to get a converse result when p 7^ 2. 

Theorem 5. For the inverse Potts problem (1) associated with 
the matrix A and data b we consider the sparsity problem 
associated with the matrix B = A'V"'' and data b' , where 
A\ V"*" and b' are given by (30), (25) and (32), respectively. 
Let u* be a minimizer of the sparsity problem with respect to 



u G arg mm 7 u 



\Bu-b'\\l 



(4) 



/i(V^u*)e (with /i given by (28)) is a 
solution of the inverse Potts problem (1) associated with A, b. 

The proofs of Theorem 4 and Theorem 5 are given in 
Subsection VI-B. 

III. Minimization of the Potts functional using 

THE alternating DIRECTION METHOD OF MULTIPLIERS 

In this section we approach the inverse Potts problem (1) 
by an iterative procedure. 

A. A new ADMM algorithm far the inverse Potts problem 

The inverse Potts problem is equivalent to the bivariate 
constrained optimization problem 

Minimize j\\S/u\\o + \\Av - b\\P 

subject to u — V = Q. 

We incorporate the constraint u — v into the target functional 
to obtain the unconstrained problem 

L^(w,u,A) =7||Vm||o + (A,M- v) 



(5) 



\Av-b\\l 



(6) 



The parameter ^ > regulates the coupling of u and v. 
The dual variable A is an n-dimensional vector of Lagrange 
multipliers. Equation (6) is called the augmented Lagrangian 
of (5). Completing the square in the second and third term of 
(6) yields 



L^{u,v,X)=-i\\Vu\\o-^\^ 
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(7) 



In order to minimize the augmented Lagrangian (7) we use the 
alternating direction method of multipliers (ADMM), see e.g. 
[10]. In the ADMM iteration we first fix v and A and minimize 
Lf^{u,v,X) with respect to u. Then we minimize L^{u,v,X) 
with respect to w, keeping u and A fixed. The third step is the 
update of the dual variable A. Thus, the alternating direction 
method of multipliers for the inverse Potts problem (1) reads 



'^ e argmin7||Vu||o + |||u - (w'" - 
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(8) 

The crucial point is that both subproblems appearing in the first 
and the second line of (8) are computationally tractable. The 
first subproblem is the minimization of a classical Potts prob- 
lem which we elaborate on in Subsection III-B. The second 
subproblem is the minimization of a classical Tikhonov-type 
problem which we explain in Subsection III-C. 

We initialize the iteration with a small positive coupling 
parameter /iq > and increase it during the iteration by a 
factor r > 1. Hence, /i is given by the geometric progression 

This assures that u and v can evolve quite independently at 
the beginning and that they are close to each other at the end 
of the iteration. We stop the iteration when the norm of m — u 
falls below some tolerance. Our approach to the inverse Potts 
problem is summed up in Algorithm I. 

Algorithm 1: iPotts-ADMM 



Input: Data 6 e R 

matrix ^GR'"''" 
Output: Computed result u G 

(1) 
begin 

V 



model parameter 7 > 0, measurement 

to the inverse Potts problem 



A"/; 
M ^ Mo; 

A^O; 
repeat 



until I 



end 



■h- Minimizer of classical L^-Potts functional (9) 
with data d = v — - and parameter S — —; 

<— Solution of Tikhonov problem (12) with data b, 
offset vector w = u + - and parameter ^ ; 

•s— A + ^{u — v); 
^— r ■ ^; 
\u-v\\'i <TOL; 



We have the following convergence result, whose proof is 
given in Subsection VI-C. 



Theorem 6. The ADMM iteration (8), and thus Algorithm 1, 
converges. 

Although we cannot expect convergence to a global min- 
imum for the NP-hard inverse Potts problem, we see in the 
experimental section that Algorithm 1 gives very satisfactory 
reconstruction results. 

In our experiments, reasonable numerical values for the 
parameters in Algorithm 1 are /io = 7-10^® as initial coupling, 
T ~ 1.05 for the increment of the coupling, and TOL == 10^^ 
for the stopping tolerance. 

B. Minimization of the classical Potts subproblem 

The first subproblem of the ADMM iteration (8) is a 
classical L^-Potts problem of the form 



P's{u)^5-\\^u\\^ + \\u-f\ 



(9) 



hi 



for parameter 5 ^ — and data / = f * 

For univariate data this problem can be solved fast and 
exactly using dynamic programming [33], [34], [35], [36]. The 
basic idea is that a minimizer of the Potts functional for data 
(/i: •••; /r) Can bc computcd in polynomial time provided that 
minimizers of the partial data (/i), (/i, /2 ),•••, (/i, •••, /r-i) 
are known. The corresponding procedure works as follows. 
We denote the respective minimizers for the partial data by 
u^, V? , ..., u^^^ . In order to compute a minimizer for data 
(/i, ..., /r), we first create a set of r minimizer candidates v^ , 
..., w*", each of length r. These minimizer candidates are given 
by 

V^ = {u'^^^,^l[l^r-],-,l^[l.T]), (10) 

^ ^ > 

Length r-l+1 

where u" is the empty vector and H[e,^r\ denotes the mean value 
of data /[£^] — {fe,...,fr). Among these candidates v^, one 
with the least Potts functional value is a minimizer for the 
data /[I,,]. 

The dynamic program for the classical Potts problem (i.e., 
the recursive computation of m" using (10)) needs at least 
O(n^) steps. The decisive step to obtain complexity 0{n^) 
is the following observation [35]. The functional value of 
a minimizer Pg{u^) for data /[i,r] can be computed from 
the functional values Pg{u^), ..., Pg{ 



of the minimizers 



u\ 



,r-l 



u' " and the mean deviation of data /[i,r]; ■••; /[r.r]- 
Indeed, using (10), the Potts functional value of the minimizer 

S)hy 



u^ is given (setting Pg{vP) 



Psi^l = mill 

£—l,....i 



P'Au' 



+ \\.f[e.,r] ~ (M[£,r], •••,Ai[£x])ll 



(11) 



The minimizing index indicates the right most jump location 
at step r and is stored as i{r). The jump locations of a 
minimizer u^ are j(r), j{j{r)), j{j{j{r))), ..., and the value 
of u^ between two consecutive jumps is given by the mean 
value of data / on this interval. Note that during the iteration 
we only have to compute and store the minimizing index j{r) 
and Pg{u^); the reconstruction only has to be done once for 
u" at the end and is thus uncritical. In total, the algorithm for 



the classical (univariate) Potts problem has 0{n) space and 
O(n^) time complexity [35]. 

For multivariate data, we cannot solve the first subproblem 
of our ADMM algorithm exactly in reasonable time because 
the classical Potts problem (9) is NP-hard in two dimensions 
[15]. However, there exist well-working practical approaches 
based on graph cuts. We here use the max-flow/min-cut based 
algorithm of the library GCOptimization 3 . [15], [37], 
[38]. 

C. Minimization of the Tikhonov subproblem 

The second subproblem of the ADMM iteration (8) is a 
classical Tikhonov problem with L^ data fitting of the form 

w\\l + \\Av - b\\P ^ mm, (12) 



A£fc I 
2 I 



where the offset vector w is given by w = ■u'^+^ + — . 

For p = 2, its solution is explicitly given by the solution of 
the normal equation 



{A*A- 



tf [d)v 



^?/+i + iA'=+A*6. 



(13) 



Here id denotes the identity matrix and A* denotes the 
transposed of the conjugate. As the time complexity of solving 
(13) is 0{n^) in general, the solution of (13) is the most 
expensive step in the ADMM iteration since the classical 
univariate Potts problem is in 0{n^). However, if A* A is 
a bandmatrix or if A* A can be diagonalized efficiently then 
the system (13) can be solved fast and we are thus able to 
deal with large data sizes. For instance, if Ax describes the 
(circular) convolution of x with some vector h, i.e.. Ax ^ h*x 
then the solution of the normal equation is given by 



\h\^ + tfj 

where r denotes the right hand side of (13). 

For p = 1 the minimization of the Tikhonov problem 
(12) is more challenging because the L^ data term is not 
differentiable. Nevertheless, the problem can be treated by 
convex optimization. We use the approach proposed in [16]. 
There, the dual problem of (12) is solved iteratively by a 
semismooth Newton method, which converges superlinearly. 
The time complexity of every iteration depends on the number 
of measurements since an m x m linear system is solved in 
each iteration. 

IV. Applications to jump-sparse recovery and 

NUMERICAL EXPERIMENTS 

In this section, we apply the inverse Potts ADMM (Algo- 
rithm 1) to the reconstruction of jump-sparse signals from 
blurred, noisy data. We consider both reconstruction from 
Fourier data and deconvolution under Gaussian, Laplacian 
or impulsive noise. (We refer to Appendix A for a formal 
description of the noise models.) We compare the results with 
the minimizers of the total variation (TV) problem given by 



j\\Vu\\, + \\Au-m 



(14) 



For the solution of this convex problem, we use the primal- 
dual method of [8] with 10 000 iterations. 



A. Deconvolution of blurred incomplete data contaminated by 
Gaussian and non-Gaussian noise 

Here, the measurement matrix A models the convolution 
with some kernel h = (h^r, ■■■,hct, ■■■,hr) of non-vanishing 
mean. We assume that only m measurements {ji, ■■■,jm}, 
m < n, are given. Hence, ^ is a reduced m x n Toeplitz 
matrix of the form 



A 



j,k 



hk-j, if \k- j\ < r 
0, else. 



(15) 



where j — ji, ...,jm, and k =- 1, ...,n. In our experiments, h 
is a Gaussian convolution kernel of standard deviation 6. 

In Figure 1, data h = Ax is corrupted by Gaussian, 
Laplacian and impulsive noise (from top to bottom) and 
TO = ^ random measurements are available. The noise 
variance is cr = 0.05 for Gaussian and Laplacian noise; in 
the impulsive noise case, 30% of the convolved signal is set 
to a random value between and 1 (uniformly distributed). For 
data contaminated by Gaussian noise we use the I? data term, 
and for the other cases the i^ data term. In the experiment 
(Figure 1) we observe that the inverse Potts ADMM algorithm 
performs as well as the total variation for impulsive noise. 
For Gaussian and Laplacian noise, the minimizers of the 
total variation problem have additional plateaus as well as 
transitional points between the plateaus. In contrast, the iPotts- 
ADMM algorithm almost perfectly recovers the jump-sparse 
signal, and, in particular, the correct number of jumps. 

B. Reconstruction of jump-sparse signal from noisy and in- 
complete Fourier spectrum 

We measure an incomplete set of m frequency components 
of a jump-sparse signal x E M". Hence, our measurement 
matrix is a reduced {m x n) Fourier matrix of the form 



-27Tijk/n 



where k = I, ...,n and j belongs to a set of m indices between 
1 and n. Such reconstruction problem have been considered 
for example in [39], [40], [41]. Here, we measure every second 
frequency component, i.e., j — 2,4, ...,n. We further assume 
that the complex valued Fourier data is corrupted by additive 
noise, i.e., 

b = Ax + T]a + 11]'^ 

where ri„,rj'^ are an m-dimensional vectors of i.i.d. Gaussian 
random variables of variance a. 

In Figure 2, we compare the performance of the inverse 
Potts algorithm (Algorithm 1) with that of TV minimization 
(14). We see that our method yields significantly higher peak 
signal-to-noise-ratios (PSNR) than minimizers of the total 
variation problem. The PSNR is given by 

PSNR(a;) ^ 10 logio ("pSfe ) ^l^) 

where x denotes the groundtruth. We further observe that 
minimizers of the total variation problem have small varia- 
tions within the plateaus and underestimate the jump heights 
("contrast reduction"). The proposed inverse Potts ADMM 
algorithm reconstructs the original signal almost perfectly. 
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Fig. 2: Reconstruction of a jump-sparse signal using only every second frequency of the Fourier spectrum. Data is corrupted by Gaussian noise (cr = 0.05). 
Tile peak signal-to-noise-ratio of inverse Potts reconstructions are significantly higher than those of the TV reconstructions. The two plots on the righthand 
side show the reconstruction results corresponding to the optimal regularization parameter with respect to the PSNR (7 = 0.21 for TV and 7 = 0.02 for 
iPotts). 
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Fig. 3: Restoration of a geometric image (dimension 256 X 256), convolved 
by a Gaussian kernel of standard deviation 12 and corrupted by Gaussian 
noise {<j = 0.05). The solution of the inverse TV functional contains smooth 
transitions between and within the geometric figures. The iPotts method 
produces no such transitions; in particular, the iPotts result is a piecewise 
constant image. (PSNR of TV: 26.829, PSNR of iPotts: 27.674). 



C. Reconstruction and segmentation of blurred images 

We use the inverse Potts functional in two-dimensions for 
the reconstruction of cartoon-Uke, i.e., piecewise constant, 
images. Such images appear for example in computed tomog- 
raphy [42]. In Figure 3, we reconstruct a cartoon-like image 
from blurred and noisy data. As in the one-dimensional setting, 
the result of the total variation minimization contains smooth 
transitions within the geometric objects and at the boundaries. 
In particular, the reconstruction is not piecewise constant. Our 
approach, in contrast, recovers the piecewise constant image 
up to rounding off the corners. 

For natural images, the Potts functional is classically used 
for (multi-label) segmentation [14], [15]. (The Potts problem 
is sometimes called the piecewise constant Mumford-Shah 




Classical Potts Segmentation 
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Fig. 4: Segmentation using the classical (9) and the inverse Potts functional (1) 
of a blurred and noisy natural image (size 241 X 161, image source [43]). Due 
to the blurring, the classical Potts segmentation exhibits additional segments at 
the boundaries, whereas segmentation with our inverse Potts ADIVIM detects 
sharp boundaries. 



problem.) We see in Figure 4 that the inverse Potts functional 
(1), which incorporates the blurring operator A, performs 
better than the classical Potts functional (9) for this task. 
Here, we segment a blurred and noisy image using the inverse 
and the classical Potts functional. Due to the blurring, the 
segmentation using the classical Potts model introduces extra 
segments at the boundaries. Minimizing the inverse Potts prob- 
lem, in contrast, detects sharp boundaries without producing 
additional boundary segments. 

V. Applications to sparse recovery and numerical 

EXPERIMENTS 

Theorem 4 asserts that solutions of the inverse Potts problem 
associated with AW yield solutions of the sparsity problem 

S^{x) = 7||a;||o -I" \\Ax - b\\P -^ min. 

Thus, we may apply the inverse Potts ADMM (Algorithm 1) 
to the sparsity problem. The corresponding method is depicted 
in Algorithm 2. 



Algorithm 2: iPotts-ADMM for the sparsity problem 

Input: Data b £ R"", model parameter 7 > 0, measurement 

matrix yl £ R""^" 
Output: Computed result x £ R" of the sparsity problem (2) 
begin 

y <— Solution of iPotts-ADMM (algorithm 1) with matrix 

AV , data /, and model parameter 7; 
X <— Vy; 
end 



We compare our method (Algorithm 2) with the following 
approaches to sparse recovery problems, which include the 
state-of-the-art methods. 

• Basis pursuit denoising (BPDN) is the convex optimiza- 
tion problem 



7l|a;||i + 



\Ax - 



b\\; 



For the experiments, we use the toolbox YALLl [31]. 
Iteratively reweighted £^ minimization (IRLl) [32] solves 
a sequence of constrained optimization problems 



|a;||i.^ 



s.t. \\Ax - 6II2 < 5, 



(17) 



where ||a;||i^^ — J2i'^i\^i\ i^ ^ weighted i^ norm. 
The weights are initialized by w^ = 1 and are updated 
depending on the solution of the previous iteration by 
"^i = jxi"- We perform five iterations and choose 
e = 10"^. We use the toolbox YALLl [31] for the 
minimization of (17). 

• Orthogonal matching pursuit (OMP) [30] greedily 
searches for minimizers of the constrained formulation 
of the L^ sparsity problem 

min ll^la; — fe||2, s.t. ||a;||o < fc. 

We use the implementation OMP . m of Stephen Becker 
available at Matlab's file exchange. 

• Iterative hard thresholding [17] uses surrogate functionals 
(forward backward splitting) for the sparsity problem. 
We here use the two variants hard_10_reg.in (IHT- 
R) and hard_10_Mterm.m (IHT-M) of the toolbox 
sparsify 0.5. 

• An ADMM method based on a "direct" splitting of (2) 
which we explain in Subsection V-B. 

A. Reconstruction of noisy and blurred sparse signals 

Our goal is to reconstruct sparse signals from noisy, blurred 
and incomplete measurements. We model this reconstruction 
task by (2) where A is a reduced Toeplitz matrix. In our 
experiments, data is blurred by a Gaussian kernel and to = ^ 
measurements are taken. Thus, we are in the setup of Sub- 
section IV-A except that now the underlying signal is sparse 
instead of jump-sparse. 

Our first example is the reconstruction of blurred and 
incomplete data under Gaussian noise (Figure 5). The noise 
distribution suggests to employ the L^ data penalty. In the 
experiment, basis pursuit denoising (BPDN) underestimates 
the height of the spikes, the Lagrangian variant of iterative 
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Fig. 6: Approximation error in dependence of the number of non-zero entries 
of solutions computed by state-of-the-art algorithms. Our iPotts-based method 
(Algorithin 2) yields lower approximation errors for any number of jumps. 
Here, the average values of 100 experiments are depicted. 



hard thresholding (IHT-R) reconstructs too many non-zero 
entries and the "direct" splitting (Subsection V-B) has to 
many additional non-zero entries. Orthogonal matching pursuit 
(OMP), iteratively reweighted i^ minimization (IRLl), hard 
thresholding (IHT-M) and the proposed iPotts-ADMM based 
approach approximate the original signal quite well; in partic- 
ular, they reconstruct the precise number of non-zero entries. 
Towards a deeper comparison of these four algorithms we 
quantify the reconstruction quality by looking at the average 
approximation error 1 1 Ax — / 1 1 2 in dependence on the number 
of non-zero entries ||x||o of a solution x; cf. Figure 6. Here, 
the average values of a series of 100 runs, using the setup of 
the experiment in Figure 5, is depicted. We observe that the 
iPotts based solutions have the least approximation errors in 
average. 

In Figure 7, we drive the same experiment as in Figure 5 
replacing Gaussian noise by impulsive noise. Due to this noise 
model, we employ the L^ data term for our iPotts-based 
algorithm. For the other methods we also use the L^ variant 
whenever it is available; to the best of our knowledge, this is 
the case for basis pursuit denoising and the direct splitting 
(Subsection V-B). We observe that the proposed algorithm 
yields an almost perfect reconstruction also in presence of 
impulsive noise and that it performs significantly better than 
the other methods in this case. 

B. Comparison with a "direct" ADMM approach to the spar- 
sity problem 

In analogy to (5), we consider the consensus form of the 
sparsity problem 



min. 



s.t. u — V = 0. 



(18) 



(19) 



7||u||o + ||At;-6||^ 

This leads to the augmented Lagrangian 

7||w||o + {\,ti-v) 
+ !^\\u-v\\l + \\Av-h\\l- 

Proceeding as in Section III we obtain a "direct" ADMM 
algorithm for the sparsity problem. This algorithm is given by 
replacing |jVu||o by ||u||o in the first line of (8). This leads 
to alternately solving a hard thresholding problem (instead of 
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Fig. 5: The original signal (dashed stem plot) is blurred by a Gaussian kernel of standard deviation 5 and corrapted by Gaussian noise of standard deviation 
CT = 0.05. We took m = ^ = 128 measurements. Orthogonal matching pursuit (OMP), iteratively reweighted £^ minimization (IRLl) and our iPotts based 
approach have the best reconstruction quality (with respect to visual inspection), followed by iterative hard thresholding (IHT-M) and the "direct" splitting (cf. 
Subsection V-B). Basis pursuit denoising (BPDN) underestimates the heights of the spikes and the Lagrangian variant of iterative hard thresholding (IHT-R) 
reconstructs too many non-zero entries. 
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Fig. 7: The same setup as in Figure 5 replacing Gaussian noise by impulsive noise (25% of data were set to a random value between —0.5 and 0.5). For 
the direct splitting (cf. Subsection V-B), basis pursuit (BPDN), and the iPotts based method we use L^ data terms. We see that the proposed L^-iPotts based 
algorithm performs significantly better than the other methods. It is able to recover the original signal almost perfectly. 



a Potts problem) and a classical TiMionov problem associated 
with matrix A. 



The difference between the "direct" ADMM approach and 
our iPotts-ADMM based method (Algorithm 2) is that they 
are based on different augmented Lagrangians. Indeed, when 
applying the iPotts-ADMM to the sparsity problem, we con- 
sider the inverse Potts problem associated with AV instead 
of A. Then, the augmented Lagrangian of the corresponding 
problem is obtained by replacing A by AW in (6). With the 



substitutions Wu = u' and Wv 



equation (6) reads 



7|k'||o + (A,V+K-i;')> 
+ l\\V+{u'-v')\\l + \\Av'-br^ 



(20) 



Comparing (19) and (20), we see that the direct method 
couples u and v directly whereas the iPotts based method 
involves the antiderivatives of u and v. 

From the experiments (Figures 5, 7, 8) we conclude that the 
iPotts-based method (Algorithm 2) is advantageous over the 
direct ADMM. In particular, the solutions of the iPotts-based 
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7 

Fig. 8: Final total energy S-y{x) of the sparsity problem using the iPotts based 
method (Algorithm 2) and the "direct" ADMM method (Subsection V-B) for 
different parameters 7. The iPotts-ADMM algorithm reaches lower energies 
over the whole parameter range. This indicates the superiority of the first 
method. Computed values are averages over 50 experiments. 



Original 




Fig. 9: Sparse image of size 30 X 30 blurred by a 7 X 7 Gaussian kernel 
(standard deviation 1) and corrupted by Gaussian noise (cr = 0.05). Data 
consists of 50 % randomly selected pixels of the middle image. The iPotts 
based method recovers almost all spikes correctly. (For visualization purposes, 
the contrast of the middle image was increased.) 



method have lower energy than the "direct" method for the 
whole range of parameters 7; cf. Figure 8. 

C. Sparse image recovery 

We also use our method to reconstruct sparse images. One 
may think of an image of small particles or of an astronomic 
image. We apply our procedure to images by reshaping the 
image to a vector and adapting the matrix A accordingly. 
Figure 9 shows the deconvolution of a sparse image using our 
iPotts-ADMM based method (Algorithm 2). In the experiment 
we see that almost all spikes are recovered while only few false 
positives ai-e reconstructed. 

VI. Proofs 

Here we provide the proofs of the theorems stated in the 
course of this paper. 

A. Existence of Minimizers 

We start out showing Theorem 1 which asserts that the 
inverse Potts problem (1) has a minimizer 

Proof of Theorem 1: In order to deal with the general case 
of a (possibly) singular matrix A we decompose the domain 
into ker A and a corresponding algebraic complement U. This 
means that U + kcrA = M" (or C") and U nVei A ^ {0}. 
(In the following we proceed without drawing attention to C" 
when writing M", but the arguments work for the complex 
case as well.) For x E M" we frequently use the decomposition 



X ^ u + V, where u — Qux is the projection Qjj of x to [/, 
and V is the corresponding projection onto ker A. 

The matrix A restricted to the subspace U is invertible, and 
since we are in finite dimensional space, there is a positive 
constant c such that 

\\Au\\ > c\\u\\ for any u G U. 

(Due to the finite dimension all norms are equivalent and 
the above inequality holds for any norm.) As a consequence, 
whenever, for a sequence u^ in U, the norm \\uk\\ tends to 00, 
the inverse Potts functional P^{uk) defined by (1) tends to 00 
as well. Therefore, for any sequence of vectors x^ in M" (not 
only in U,) we obtain the implication: 

P^{xk) is bounded 

=^ Uk = QuXk has a converging subsequence. (21) 

This is a consequence of Ax}. — AQ(jXk- 

Our next preparatory step introduces the mapping s on U 
which assigns to each u E U the minimal number of jumps 
of all vectors in u + ker A, i.e.. 



s{u) 



mm 



\Vv\ 



We show that this mapping s is lower semicontinuous which, 
in our context, means that the preimages of the sets {0, . . . , fc} 
are closed for all /c € N. To see this, we first observe that the 
set Affc of all vectors in M" with at most k jumps is structurally 
a finite union of vector spaces (of dimension k + 1.) More 
precisely, 

Mk - {x e M" : ||Vx||o < fc} = U Xj, 

JC{l,...,n-l},|J| = fe 

where Xj are those vectors whose jump sets are contained 
in Jc {!,...,«— 1}. Furthermore, a vector u G U has the 
property s{u) < fc if and only if there is a vector x E Mk 
(i.e., with at most fc jumps) such that Qi/x = u. Summing up, 

5-i({0,...,fc}) = Q[/(Affc)= y QuiXj). 

JC{l,...,n-l}.|J|=fc 

(22) 
We discuss the right hand side of (22) to see the lower 
semicontinuity of s. Each Qu{Vj) is a finite dimensional linear 
subspace and thus closed; so as a finite union of closed sets 
the right hand side of (22) is closed. Therefore the left hand 
side of (22) is closed which by definition implies the lower 
semicontinuity of s. 

Now we can show the assertion of the theorem. We consider 
a sequence Xk such that the values P^{xk) of the inverse Potts 
functional P^ tend to an infimum, i.e., 

lim P^{xk) = inf P^(x). 

For every member of the sequence, we write Xk = Uk + Vk 
with Uk E U and Vk E ker A. By (21) we find a subsequence 
Xki such that Uki — Qu^ki converges to some u E U. Since 
P^{xki) converges and Auki ~ Axk^ we have that 

|||Va;fc,||o-||Va;fcJ|o| < \P.,{xk,) ~ P-yixkJ\ 

+ \\\Auk, -b\\P- \\Auk,. -b\\^p\ ^ as /,r^- 00. 
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This means that, for sufficiently large I, the number of jumps 
llVxfe, llo becomes constant; let us denote this constant by j. 
As a consequence, s(uj.J < ||VxfeJ|o — j, and thus, by the 
lower semicontinuity of s, s{u) < j. Hence, by the defintion 
of s, there is a vector x* E u + kci A such that the number 
of jumps of X* is smaller than or equal to j. Then, 



P^{x* 



Vx* 



\Ax* 



b\\; 



Vx* 



lAu- 



b\\l 



< 



j + \iiR\\Auki -b\\P = limP^(a;A;,) 



which shows that x* is a minimizer as desired. ■ 

Next we show Theorem 2 which states that the continous- 
time analogue of Theorem 1 is wrong. We give counterexam- 
ples, i.e., we find bounded operators A and data / such that 
the continuous-time inverse Potts functional 

P,{u)=j-\\\/u\\o + \\Au-f\\;, 

if u is a piecewise constant function on the interval [0, 1], and 
P-y{u) ~ oo Otherwise, has no minimizer 
Proof of Theorem 2: 
We consider a positive function g E LP[0, 1] with total mass 
1 which is supported in the interval [| — e, | +&] with positive 



Our 



£ < |. We use the symbol g for its left-shift by i. 
counterexamples are the (cyclic) convolution operators with 
functions g as above, i.e., operators A defined by Au = g* u, 
and the data given by f = g. 

We claim that, for Potts parameter 7 with 7 < 70 (defined 
in (24) below,) 



ini P^{v) ^ 2j < P^{u) for all u. 



(23) 



This means that there is no minimizer in that case and thus 
shows the assertion of the theorem. In order to show the 
equality in (23), we consider the sequence of characteristic 
functions ?i„ = f Ifi-i - + -]■ ^^ have that ||Vu„||o ~ 2 
and \\Aun — f\\ -^ 0. Thus, P^{un) -^ 27. This yields 
inf^,P^(f) < 27. It remains to show the inequality in (23) 
(which in turn implies the equality in (23).) To this end, we 
have to consider the set of functions u with at most one jump 
and find 7 > such that d{u) = || Au- /||^ > 27 for all such 
u. The set B = {x E [0, 1] : f{x) > 2} has positive Lebesgue 
measure X{B) since / has total mass 1 and is supported on an 
interval of length bounded by j.If u < 1, then d{u) > MB). 
So in order to obtain d{u) < \{B), we need that m > 1 
either to the left or to the right of its (sole) jump location. 
Then Au > 1 on at least one of the intervals [e, | — e] and 
[| + £, 1 — e]. Both of these intervals have length ^ — 2e, 
and, on both intervals, / = 0. Therefore, if d(u) < X{B), we 
necessarily have d{u) > ^ — 2e > ^ . Then, for any u with at 
most one jump. 



d(M)>min(A(B),i)=:27o. 



(24) 



If u has two or more jumps then trivially P^{u) > 2j. 
Together, this implies that, for any 7 with 7 < 70, the inverse 
Potts functional P^ fulfills P^{u) > 27 for all u G L'p[0, 1]. 
This shows (23) which completes the proof. ■ 



B. Relations to sparsity 

We first show Theorem 4 which shows how to transform a 
sparsity problem into a jump-sparsity problem. 

Proof of Theorem 4: For x* satisfying (3), we define 
u* = Vx*. Towards a contradiction we assume that there is 

ueK" such that 7||u||o + Pu-6||P < j\\u*\\o + \\Au* -b\\P. 
Then, for u, there is a; G M"+^ such that u = Vx. Then, 

7||Vx||o + \\AVx - b\\P < j\\u*\\o + \\Au* - b\\P 

^^\\Vx*\\o + \\AVx*-b\\P. 

which is a contradiction. ■ 

For p = 2 we show a converse statement. It is formulated 

as Theorem 5 and proved next. In its proof we make use of 

the decomposition of M" into the orthogonal direct sum M" = 

V Me, where e denoted the constant vector (1, . . . , 1)-^ and 

V is the linear space of vectors with zero mean. Observing 
that the linear operator V is bijective from the linear space V 
to K"^^, we use the symbol V+ for the mapping M"^^ -^ V, 



I" 1, we use the symbol V+ for the mapping M" ^ 



(25) 



for the inverse of the mapping V restricted to the subspace V. 
Proof of Theorem 5: We consider the inverse Potts 
functional given by (1) for p = 2. We decompose x E M" 
according to x = Xq+x, with xq E V,x G Me. Applying this 
decomposition to (1) yields 

P-yi^) = 7l|V(xo + x)||o + II Axo + Ax-b\\l 



7||Vxo||o + ma;o + Ax-6|| 



2- 



We write x = /ie to obtain 



P^(x) ^ 7l|Vxo||o + ||Axo + nAe - b\\ 



(26) 



Let us fix xq for the moment and let us look for /i = ^(^^o) 
which minimizes the function /i -^ P-y{xQ + fie). Since 
||V(xo + Aie)||o = ||V(xo + M'e)||o for all fi,fi' we have to 
minimize (w.r.t. /i) 



^ I /x ^ Ay + ^ A^jxo.j - bi 



(27) 



The corresponding minimizer /i(xo) can be computed explic- 
itly (e.g., by derivating). It is given by 

Ai(a;o) = ^2 ' (28) 

where Ai is the sum of the i*'* row of the matrix A given by 
(31). In particular, /i(xo) depends affine linearly on xo, i.e. is 
of the form d — Exq where d is a constant and _E is a row 
vector of length n, both not depending on xq. Plugging the 
expression (28) for /i(xo) into (26), we obtain a minimization 
problem in xq. It is given by 

r,'l|2 , „• 



7l|Vxo||o + ||A'xo 



where A' is the matrix given by 



Akj — A-kj 



YZX 



(29) 



(30) 



with 



Ai :— 2_^ Aij 



and b' is the vector given by 



61 






(31) 



(32) 



After these preparations we show the theorem; we consider 
a minimizer u* of the sparsity problem (4) w.r.t. the matrix 
B ~ A'y and data b' . The crucial point is that V is an iso- 
morphism from V onto M"^^ which implies the equivalence 

u* minimizes (4) <^ Xq = V^m* minimizes (29) . 

Applying this equivalence, Xq = V+u* is a minimizer of (29), 
and, using (28), the vector x* = Xq + l^i^'o) is a minimizer 
of the original Potts problem (1) for A,b. ■ 

Using the relation between inverse Potts and sparsity prob- 
lems we are now able to show the complexity statement 
Theorem 3 which asserts NP-hardness of the inverse Potts 
problem. 

Proof of Theorem 3: The sparsity problem (2) is NP- 
hard by [26, Theorem 3] (p > 1, 7 > 0.) According to 
Theorem 4 each instance of the sparsity problem (2) defines an 
instance of the inverse Potts problem (1). In particular, for any 
NP-hard instance of the sparsity problem (with matrix A and 
data b) there is a corresponding inverse Potts problem (with 
matrix AW and data b.) The transformation of the functionals 
and the transformation of the corresponding minimizers given 
by Theorem 4 can obviously be done in polynomial time. 
Therefore the Potts problem is NP-hard. ■ 

C. Convergence 

In our presentation we have assumed that the sequence /i^ 
is a geometric progression. What we actually need is that /ij. 
is a non-decreasing sequence fulfilling 



V^= < 00 



(33) 



which is obviously satisfied for geometric progressions. So we 
show Theorem 6 assuming (33) instead. 

Proof of Theorem 6: We consider the Potts ADMM 
iteration for u'^, v'^ and A'^ given by (8). We show that 



{u^ ,v'') -^ {u* ,v*) with u* = V* 



and 



0, (34) 



consequence of the minimality of m''+^. Using the definition 
of Fk yields 



,fe+l|| , Mfen fc+l 



7l|V«^+^||o+^||«^+^-h -- I 



A* 



Mfe 



A'= 



<7||VU'= ||o<7", 

V A*fc/ 

where n is the length of Vk- Since the first summand on the 
left hand side is non-negative we get that 



In particular. 






lim u^+^ -[v'' - — \=Q. 

^J■k, 



(35) 



(36) 



Now we draw our attention to the second line of (8). We 
denote the corresponding functional by 

G.W = p. -6||^+f||. -(,/+! + £) 11^. 

The minimality of v^^^ implies Gk{v^^^) 
< Gk (u''+^ + ^) • We apply the definition of Gk 
and estimate 

\\Av''+^-b\\P+^\\v'+'- \u 



Hk 



<\\A u 



,k+l 



A* 






br. 



<\\A(u''+^ + — - v'') + Av'^ - b\\P 



<(\\A\\\\u''+' + — -v''h + \\Av''-b\\p 

(37) 

Here ||A|| is the norm of A viewed as an operator from £^ 
to iP. We combine the inequalities (37) and (35) in order to 
obtain that 

\\Av''+'-b\\,<^^^^ + \\Av'-b\\,. 



^^k 



Solving this recursion yields 



which is a qualitative version of the assertion of the theorem. 
We denote the functional occurring in the first line of (8) 

by Fk, i.e.. 



^fe(u)=7||V«||o + yh 



^J■k 



Using this notation, the first line of (8) reads m'^+^ e 

argmiiiti Ffe(M). In order to estimate ||7i'^+^ — [v'^ — ~)l|2 

which is a 



\\Av^+'-b\\,<\\AhnJ2- + Uv'-b\\,, 
3=1 ^' 

which shows that the sequence (|| Af '^+^ — &||p)fegN is bounded. 
Together with (37) this implies 

\ ^J'k J ^J■k 

where C is a positive constant independent of k. Using (36) 
we get that 



observe that Fk{u''+^) < Fk 



,k+l _ 



,,fe+l I „,k + l , ^ \ ||2 



Mfe 



fik\\v - [u 



is bounded. (38) 
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We show the convergence of the sequence v^ by showing 
that it is a Cauchy sequence. To this end we estimate 



k'^+i-v'^l 



< Ib'^+i - u'^+i - 



-I 



,k+l 



A* 
Mfe 



-v% 



Now we apply (35) and (38) which yield 

C 



,,k+l 



< 



for some constant C > which is independent of k. Assump- 
tion (33) on /ife guarantees that v'^ is a Cauchy sequence and 
hence that v'^ converges to some v* . 

We use the third line of (8) to obtain the equality 



Mfc+i Mfe+i 



A'^ 



(— +w'=+^-i;'=) + («'=- 1;''+^) 



fe+i- 



Mfe 



(39) 



By (36) and (38) each term in parenthesis converges to 0. 
Since /i^ is non-decreasing, we have that ^k/ i^k+i < 1 and, 
therefore, (39) implies that 



lim — 



0. 



We rewrite the third line of (8) as m'^+^ 



„fe+i 



(A' 



fe+i 



\^)/ ^ik and use the monotonicity of the sequence fXk to obtain 
the inequality 



,fe+i 



,,fe+i| 



< 



IIA''^ 



A''^ 



< 



A''^ 



Mfc+i 



This means that u^ — v^ ^^ Q and, since v^ converges, also 
Uk converges and the corresponding u* equals v* . This shows 
(34) and completes the proof. ■ 

VII. Conclusion and outlook 

We have shown that the inverse Potts problem has a 
minimizer in the discrete setting but that the time continuous 
counterpart does not have minimizers in general. We further 
have shown that the computation of minimizers is an NP- 
hard problem. Having accepted that the computation of exact 
solutions are unfeasible, we have proposed a new approach to 
the inverse Potts problem based on the alternating direction 
method of multipliers. In our experiments we have compared 
the iPotts-ADMM algorithm with total variation minimization 
for jump-sparse reconstruction. We have observed that our 
method often performs better than but at least as well as 
TV minimization. We further have shown that the sparsity 
problem can be reduced to an inverse Potts problem for p > 1 . 
The experiments indicate that the iPotts-based approach to the 
sparsity problem performs as least as well as the state-of-the- 
art algorithms in presence of Gaussian noise and significantly 
better in presence of impulsive noise. 

Future research aims at faster algorithms for the multivariate 
inverse Potts problem and at Potts problems with manifold 
valued data. 



Appendix A 

We consider Gaussian, Laplacian, and impulsive noise. The 
first two types of noise are additive. Thus the measurement is 
given by 

b = Ax + rj„, 

where 77^ is a ?7i-dimensional vector of i.i.d. random variables 
of standard deviation a. In case of Gaussian noise, the prob- 
ability density function is given by 

1 _ x2 

In the case of Laplacian noise, the density is defined by 

1 _yii 



■p{x) = 



aV2 



-^N 



In the case of impulsive noise, we randomly choose a pre- 
scribed percentage of indices / between 1 and n and set each 
data point belonging to that index set to a random number, i.e. 




Here, ^ is a random variable which is uniformly distributed in 
the interval [0, 1] for the jump-sparsity experiments and in the 
interval [— |, i] for the sparsity experiments. 
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